! This file is part of PanGu. All rights reserved.
module pangu

    use seakeeping_kinds, only: rk
    use seakeeping_math, only: arange
    use seakeeping_constants, only: fmt
    use seakeeping_filesystem, only: operator(.join.), is_exist, mkdir
    use pangu_D_Alembert, only: get_acceleration, D_Alembert, get_center, get_moment_of_inertia
    use pangu_write, only: write_load_file, write_acc_file
    use pangu_read, only: read_pangu_file, read_data_file

contains

    !> 计算剖面载荷
    subroutine get_section_load(pangu_file, output_dir, load_file, acc_file, error)
        character(*), intent(in) :: pangu_file, load_file, acc_file  !! 配置文件、质量文件、力文件、载荷文件、加速度文件
        character(*), intent(in) :: output_dir  !! 输出目录
        character(:), allocatable, intent(out) :: error  !! 错误信息
        integer :: nstation  !! 站数
        real(rk), allocatable :: mass(:), pressure(:)  !! 质量、外力
        real(rk) :: acceleration(2)  !! 加速度
        real(rk) :: iy, xg, width  !! 纵摇转动惯量，重心位置，船宽
        real(rk), allocatable :: load(:, :), acc(:)  !! 载荷：剪力、弯矩，加速度
        real(rk), allocatable :: x(:)  !! 站中心位置
        real(rk), parameter :: dx = 1.0_rk  !! 站距

        call read_pangu_file(pangu_file, nstation, width, error)
        if (allocated(error)) return

        allocate (mass(nstation), pressure(nstation), load(2, 0:nstation), acc(nstation))

        call read_data_file(pangu_file, mass, nstation, 'mass', error)
        if (allocated(error)) return
        mass = mass*1000.0_rk  ! 质量单位转换：t -> kg

        call read_data_file(pangu_file, pressure, nstation, 'pressure', error)
        if (allocated(error)) return

        allocate (x, source=arange(0.5_rk, (nstation - 1) + 0.51_rk, 1.0_rk)*dx)
        call get_center(mass, x, xg)
        print fmt(1), 'Xg(m): ', xg

        call get_moment_of_inertia(mass, x, xg, dx, iy)
        print fmt(3), 'Iy(N*m): ', iy, 'radius of inertia(m): ', sqrt(iy/sum(mass))

        call get_acceleration(mass, pressure, xg, x, dx, iy, width, nstation, acceleration)
        print fmt(1), 'Acceleration(m/s^2, rad/s^2): ', acceleration

        call D_Alembert(mass, pressure, acceleration, xg, x, dx, width, load, acc, nstation)

        if (.not. is_exist(output_dir)) call mkdir(output_dir)

        call write_load_file(output_dir.join.load_file, load, nstation, error)
        if (allocated(error)) return

        call write_acc_file(output_dir.join.acc_file, acc, nstation, error)
        if (allocated(error)) return

    end subroutine get_section_load

end module pangu
